#include <qmath.h>

#include "qlsfft.h"
#include "fftw3.h"

QLsFFT::QLsFFT()
{

}

/* 傅里叶变换 */
void QLsFFT::FFT(const double *in, double *out, int count)
{
	if(!in) return;

	/* 向高位转换为2的倍数 */
	count = TransTwoPowerHigh(count);

	fftw_complex *fftw_in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * count);
	fftw_complex *fftw_out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * count);
	fftw_plan fftw_p = fftw_plan_dft_1d(count, fftw_in, fftw_out, FFTW_FORWARD, FFTW_ESTIMATE);

	/* 存入数据 */
	for (int i = 0; i < count; i++)
	{
		fftw_in[i][0] = in[i];		/* real */
		fftw_in[i][1] = 0;			/* image */
	}

	/* 傅里叶变换 */
	fftw_execute(fftw_p);			/* repeat as needed */

	/* 转换为频域数据的个数为时域数据个数的一半 */
	count >>= 1;

	/* 频域转换 */
	for(int i = 0; i < count; i++)
	{
		double val = sqrt(fftw_out[i][0] * fftw_out[i][0] + fftw_out[i][1] * fftw_out[i][1]);
		out[i] = 20 * log10(val);
	}

	fftw_destroy_plan(fftw_p);
	fftw_free(fftw_in); fftw_free(fftw_out);
}
void QLsFFT::FFT(const QVector<double> &in, QVector<double> &out)
{
	/* 根据输入数据个数，设置输出长度 */
	int count = TransTwoPowerHigh(in.size());
	out.resize(count >> 1);

	FFT(in.data(), out.data(), count);
}

/* 向高位转换为2的倍数 */
inline int QLsFFT::TransTwoPowerHigh(int num)
{
	num--;
	num |= num >> 1;
	num |= num >> 2;
	num |= num >> 4;
	num |= num >> 8;
	num |= num >> 16;
	num++;

	return num;
}
